
# # # # # # # # # # # # # # # # # # # # # 
#
# Appendix D: Constitutional Integration Indices
# 
# Part of replication of:
# Traditional Authorities, Norm Collisions, and Communal Conflict 
#
# Journal of Conflict Resolution
#
# By Clara Neupert-Wentz, 2024
#
# # # # # # # # # # # # # # # # # # # # # 

#Note: main replication data (nc_rep) contains all relevant variables 
#The following code shows how the constitution indices are coded (Appendix D)


###### Set-up ----

#Libraries
library(foreign)
library(tidyverse)
library(haven)
library(poLCA)
library(xtable)

#Cleanup
rm(list = ls())

#Globals
set.seed(9999)
setwd("/Neupert-Wentz Replication/")

# Load country level const_rep data
load("data/const_rep.RData")



###### Principal Component Analysis ----


# Select variables for PCA
vars_for_pca <- c("leader", "bodyack", "moreTPI", "disres_cst", "law_cst",
                  "bodyest", "off", "sanctno", "sanctyes", "negrep", "negparty", "funcnrl", "funccc", 
                  "funcint", "funcrep", "funcac", "funcas", "funco", "funcdma", "funcdmc", "funcv", "funcap", 
                  "funcjc", "funcjs", "funcdev", "funcec", "lawcoll_cst", "lawsec_cst", "conscstrel",
                  "conshril", "finstat", "fintax", "finrev", "proh", "conres", "conscon", "conssec")

df_pca <- const_rep[, vars_for_pca]

# Compute PCA
pca_res <- prcomp(const_rep[, vars_for_pca])

# Prepare the dataframe for the PCA scores
const_rep$pca_score <- NA

# Predict the first principal component scores
const_rep$pc1 <- predict(pca_res)[, 1]

# Flip sign if the loading of the first variable on the first principal component is negative
if (pca_res$rotation[1,1] < 0) {
  const_rep$pc1 <- -1 * const_rep$pc1
}

# Rescale to 0 - 1
min_val <- min(const_rep$pc1, na.rm = TRUE)
max_val <- max(const_rep$pc1, na.rm = TRUE)
const_rep$pc1 <- (const_rep$pc1 - min_val) / (max_val - min_val)

# > Print Table A2 ----
pca_result <- print(pca_res)
result_df <- data.frame(pca_result$rotation[, 1:5])
xtable(result_df)

pca_summary <- summary(pca_res)
summary_df <- data.frame(pca_summary$importance[, 1:5])
xtable(summary_df)






###### Latent Class Analysis ----


#Prepare: Recode dummies to 1/2
tpi_lca <- const_rep %>%
  mutate(across(c(leader, bodyack, moreTPI, disres_cst, law_cst, 
                  bodyest, off, sanctno, sanctyes, negrep, negparty, funcnrl, funccc, 
                  funcint, funcrep, funcac, funcas, funco, funcdma, funcdmc, funcv, funcap, 
                  funcjc, funcjs, funcdev, funcec, lawcoll_cst, lawsec_cst, conscstrel, 
                  conshril, finstat, fintax, finrev, proh, conres, conscon, conssec),
                ~{.x + 1})) 

f <- cbind(leader, bodyack, moreTPI, disres_cst, law_cst,
           bodyest, off, sanctno, sanctyes, negrep, negparty, funcnrl, funccc,
           funcint, funcrep, funcac, funcas, funco, funcdma, funcdmc, funcv, funcap, 
           funcjc, funcjs, funcdev, funcec, lawcoll_cst, lawsec_cst, conscstrel, 
           conshril, finstat, fintax, finrev, proh, conres, conscon, conssec)~1

#Run LCA for different no. of classes
one_class <- poLCA(f,tpi_lca,nclass=1) 
two_classes <- poLCA(f,tpi_lca,nclass=2)
three_classes <- poLCA(f,tpi_lca,nclass=3,maxiter=8000) 
four_classes <- poLCA(f,tpi_lca,nclass=4,maxiter=8000) 
five_classes <- poLCA(f,tpi_lca,nclass=5,maxiter=8000) 
six_classes <- poLCA(f,tpi_lca,nclass=6,maxiter=8000) 

# > Figure A1: Compare BIC ----
plot(c(one_class$bic, two_classes$bic,
       three_classes$bic, four_classes$bic,
       five_classes$bic, six_classes$bic),
     ylab = "BIC", xlab = "Number of Classes", main = "BIC for Different Number of Classes")


# > Figure A2: LCA for two classes ----
plot(two_classes, mar = c(6, 6, 6, 12), las = 1)

# Extracting latent class probabilities
tpi_lca$two_class <- two_classes$predclass

# Recode back to 0/1 (0 = little integration (154 countries) and 1 = a lot of integration (39 countries))
tpi_lca$two_class01 <- ifelse(tpi_lca$two_class == 2, 1, 0)
table(tpi_lca$two_class01, tpi_lca$two_class )


#### NA for years before last constitutional change ----

#Load main replication data
load("data/nc_rep.RData")


# Join datasets

#Note: replication dataset already contains all variables, this is for demonstration purposes only
merged_data <- left_join(nc_rep, tpi_lca, by = "cow")

## NA for years prior to last constitutional change
merged_data <- merged_data %>%
  mutate(across(c(two_class01.y, pc1.y), ~ ifelse(year < mcv_precise.y, NA, .), .names = "{.col}_mcv"))

# Check if variables are the same
table(merged_data$two_class01_mcv, merged_data$two_class01.y_mcv, useNA = "ifany")
cor(merged_data$pc1_mcv, merged_data$pc1.y_mcv, use="complete.obs")



